function bci = bciInit(bci)

global MEMepoch;
global MEMdata;

% sample sizes
bci.init.nOffset = round(bci.param.winoffset*bci.srate);
bci.init.nWinLength = round(bci.param.winsize*bci.srate);
bci.init.nBaseline = round(bci.param.baselineLength*bci.srate);
bci.init.nEpochLength = bci.init.nWinLength + 2*bci.init.nOffset; % length of Epoch in samples
bci.init.nWinShift = round(bci.srate/bci.param.winsrate);
bci.init.intervalOfInterest = bci.init.nOffset+1:bci.init.nOffset+bci.init.nWinLength;

if bci.param.specMethod == 2,
    % taper params
    bci.init.taperParams.tapers = bci.param.tapers;
    bci.init.taperParams.Fs=bci.srate;
end

if bci.param.specMethod >3,
    % exact FT parameters (fouriertransform at specified frequencies only in
    % whole period intervals)
    bci.init.eFT.periodMaxIdx = zeros(1,length(bci.param.freq));
    if any(bci.param.specMethod==[5 9 10 11]),
        nTap = max(bci.param.tapers(:,2));
    else
        nTap = 1;
    end
    freqTap(1:length(bci.param.freq)) = bci.param.tapers(:,2);
    % assume bci.nChan channels and preallocate
    bci.init.eFT.eFunc = zeros(bci.nChan,bci.init.nWinLength,length(bci.param.freq),nTap);
    eFTtapers = zeros(size(bci.init.eFT.eFunc));
    for fi = 1:length(bci.param.freq),    
        taperParam = bci.param.tapers(min(size(bci.param.tapers,1),fi),:);
        nPeriods = floor(bci.init.nWinLength/bci.srate *bci.param.freq(fi)); % number of periods
        nSamp = round(nPeriods/bci.param.freq(fi)*bci.srate);% number of sample points
        t=linspace(0,bci.param.winsize,bci.init.nWinLength);
        bci.init.eFT.periodMaxIdx(fi) = nSamp;
        eFunc = exp(-1i*2*pi*bci.param.freq(fi)*t);
        bci.init.eFT.eFunc(:,:,fi,:) = repmat( repmat( eFunc,bci.nChan,1),[1,1,1,nTap]);
        bci.init.eFT.eFunc(:,nSamp+1:end,fi,:)=0;
        tapers = permute(dpsschk(taperParam,nSamp,bci.srate),[3 1 2]);
        eFTtapers(:,1:nSamp,fi,1:freqTap(fi))=permute(repmat( tapers,bci.nChan,1), [1 2 4 3]);     
    end
    if any(bci.param.specMethod==[5,9,10]),
        % multiply taper windows with e function
        bci.init.eFT.eFunc = bci.init.eFT.eFunc.*eFTtapers;
    end       
end

if bci.param.specMethod==2,
    [dummy,centerFreq]=mtspectrumc(rand(1,length(bci.init.intervalOfInterest))',bci.init.taperParams);
    TW = bci.param.tapers(1)/bci.param.winsize;
    freqResolution = TW/3;
elseif bci.param.specMethod==0 && ~isempty(bci.param.resampFreq),
    centerFreq = bci.param.resampFreq/2*linspace(0,1,2^nextpow2(...
                    ceil(bci.init.nWinLength*floor(bci.param.resampFreq)/floor(bci.srate)))/2);
else
    centerFreq = bci.srate/2*linspace(0,1,2^nextpow2(bci.init.nWinLength)/2);
    freqResolution = centerFreq(2);
end

% determine default frequencies
if isempty(bci.param.freq),
    maxFreq=find(centerFreq>100,1);
    if isempty(maxFreq),
        maxFreq=length(centerFreq);
    end
    bci.param.freq=intersect(unique(round(centerFreq(2:maxFreq))),round(0:freqResolution:maxFreq));
end

% determine number of features
if bci.param.specMethod==0,
    if ~isempty(bci.param.resampFreq),
        bci.init.nFeatPerChan = ceil(bci.init.nWinLength*floor(bci.param.resampFreq)/floor(bci.srate));
        bci.init.resampIntOfInterest = (1:bci.init.nFeatPerChan) + ...
            ceil(bci.init.nOffset*floor(bci.param.resampFreq)/floor(bci.srate));
        bci.init.filtSrate = bci.param.resampFreq;
    else
        bci.init.nFeatPerChan = length(bci.init.intervalOfInterest);
        bci.init.filtSrate = bci.param.resampFreq;
        if isempty(bci.init.filtSrate),
            bci.init.filtSrate = bci.srate;
        end
    end
elseif bci.param.specMethod>0 && bci.param.specMethod<3 || bci.param.specMethod>3,
    bci.init.nFeatPerChan = length(bci.param.freq);
elseif bci.param.specMethod==3,
    bci.init.nFeatPerChan = bci.param.ARorder;
end

%default epoch rate:
bci.init.epochRate = bci.param.winsrate;

if bci.paradigm.id==12||bci.paradigm.id==14, % P300 case
    if bci.paradigm.doAverageFlashEvents,
        bci.init.nFeatPerChan = 2*round(bci.param.resampFreq*bci.paradigm.flashWin)...
                                * length(bci.eventsToClassify) ; % 2classes by nSamp by nTargets
    else
        bci.init.nFeatPerChan = bci.paradigm.nTargetFlashes*round(bci.param.resampFreq*bci.paradigm.flashWin)...
                                * length(bci.eventsToClassify) ; % nTargetFlashes by nSamp by nTargets
    end
        bci.init.epochRate = 1/bci.movDur;                            
end

% combined feature spaces:
if any(bci.param.specMethod == [8,9]),
    bci.init.nFeatPerChan = bci.init.nFeatPerChan*2;
end

bci.init.featOfInterest = 1:bci.init.nFeatPerChan;

% filter params
if length(bci.param.bandFreqs)==2,
    Wn_high = bci.param.bandFreqs(1)/(bci.srate/2);             % highpass Cutoff Frequenz 
    Wn_low = bci.param.bandFreqs(2)/(bci.srate/2);              % lowpass Cutoff Frequenz 
    [bci.init.filtCoeffB,bci.init.filtCoeffA] = ...
        butter(3,[Wn_high Wn_low]);                 % bandpass filter coefficients
end

% FourierParams
if bci.param.specMethod==0 && ~isempty(bci.param.resampFreq), % if time series, use resampled
    bci.init.NFFT = 2^nextpow2(bci.init.nFeatPerChan);
%    bci.init.fftFreqIdx=round(bci.init.NFFT.*bci.param.freq/bci.param.resampFreq+1);
    bci.init.fftWin = window(eval(['@' bci.param.fftWin]),bci.init.nFeatPerChan)';
else    
    bci.init.NFFT = 2^nextpow2(bci.init.nWinLength);
%    bci.init.fftFreqIdx=round(bci.init.NFFT.*bci.param.freq/bci.srate+1);
    bci.init.fftWin = window(eval(['@' bci.param.fftWin]),bci.init.nWinLength)';
end

% determine index mapping to centerFrequencies
bci.init.fftFreqIdx = zeros(size(bci.param.freq));
for k=1:length(bci.init.fftFreqIdx),
    [dummy, bci.init.fftFreqIdx(k)]=min(abs(centerFreq-bci.param.freq(k))); % find closest index
end
if length(unique(bci.init.fftFreqIdx))~=length(bci.init.fftFreqIdx),
    warning('Multiple frequencies point to center frequency vector (redundant features)');
end

% classifier function handles
if strcmpi(bci.param.classifier,'mimo'),
    bci.param.trainFunc=str2func('bciTrainMIMO');
    bci.param.testFunc=str2func('bciTestMIMO');
elseif strcmpi(bci.param.classifier,'svmmulti'),
    bci.param.trainFunc=str2func('bciTrainSVMMulti');
    bci.param.testFunc=str2func('bciTestSVMMulti');
elseif strcmpi(bci.param.classifier,'logreg'),
    bci.param.trainFunc=str2func('bciTrainLogReg');
    bci.param.testFunc=str2func('bciTestLogReg');
elseif strcmpi(bci.param.classifier,'bayeslda'),
    bci.param.trainFunc=str2func('bciTrainBayesLDA');
    bci.param.testFunc=str2func('bciTestBayesLDA');
elseif strcmpi(bci.param.classifier,'rda'),
    bci.param.trainFunc=str2func('bciTrainRDA');
    bci.param.testFunc=str2func('bciTestRDA');
else
    bci.param.trainFunc=str2func('bciTrainSVM');
    bci.param.testFunc=str2func('bciTestSVM');
end
if bci.paradigm.id==12||bci.paradigm.id==14,
    bci.param.trainFunc=str2func('bciTrainOddball');
    bci.param.testFunc=str2func('bciTestOddball');
end
% enable fast recorder system check 0= ft buffer, 1=tdt,2=btini
bci.init.recordingSystemID = double(strcmpi(bci.recordingSystem,'tdt'))+double(strcmpi(bci.recordingSystem,'btini')*2);

% default setting for noise cancellation
bci.init.doNoiseCancellation = false;

% MEG transfer information
if isfield(bci,'configFile') && ~isempty(bci.configFile)
    load(bci.configFile);%BTIconfig = HDR;
    bci.init.trigChan=HDR.trigChan;
    bci.init.respChan=HDR.respChan;
    % BTI listener on SUN machine already sorts the channels using *.cfg
    % file
    bci.init.dataChan=HDR.MEGtransferChan;
    [sortedIdx transIdx] = sort(HDR.sendIdx);
    transIdx=transIdx(sortedIdx>0);
    bci.init.channelGain=HDR.ChannelGain(transIdx);
    bci.init.channelUnitsPerBit=HDR.ChannelUnitsPerBit(transIdx);
    bci.init.calib = single(diag(bci.init.channelGain .*bci.init.channelUnitsPerBit));
    if isfield(bci,'referenceChan') && ~isempty(bci.referenceChan),
        bci.init.doNoiseCancellation = true;
    end
else
    if isfield(bci.param,'trigChan') && ~isempty(bci.param.trigChan),
        bci.init.trigChan = bci.param.trigChan;
    else %
        bci.init.trigChan=[];
    end    
    if isfield(bci.param,'respChan') && ~isempty(bci.param.respChan),
        bci.init.respChan = bci.param.respChan;
    else %
        bci.init.trigChan=[];
    end
    bci.init.calib = eye(bci.nChan);
    bci.init.dataChan =bci.preprocChan;
end

% spatial smoothing params
if bci.param.spatialFilter>0,
    chanIdx = bciFind(bci.init.dataChan,bci.preprocChan);
    bci.init.spatialFilter = bciGetSpatialFilterCoeffs(bci.param.spatialFilter,bci.param.chanDistThresh,chanIdx,bci.param.spatialMethod)';
end


% triggerport initialization
if bci.param.doSendTrigger %&& ~isfield(bci.init,'iOObj'),
    if isempty(bci.param.trigAddress), % DIO 
        bci.init.iOObj = digitalio('mcc',0);
        addline(bci.init.iOObj,0:7,'out');
    else % parallel port
        %create IO32 interface object
        if strcmp(mexext,'mexw64'),
            bci.init.iOFunc = str2func('io64');
        else
            bci.init.iOFunc = str2func('io32');
        end
        bci.init.iOObj = bci.init.iOFunc();
        % %install the inpout32.dll driver. Status = 0 if installation successful
        iOStatus = bci.init.iOFunc(bci.init.iOObj);
        if(iOStatus ~= 0)
            disp('parallel port inp/outp installation failed!');
            errordlg('Parallel port inp/outp installation failed!','Trigger port Error');
            bci.init.iOObj=-1;
        else
            bci.init.iOFunc(bci.init.iOObj,bci.param.trigAddress,0);
        end
        
    end
end

% load sound file
if isfield(bci.paradigm,'fbSoundFile') && ~isempty(bci.paradigm.fbSoundFile),
    load(bci.paradigm.fbSoundFile);
    bci.init.feedbackSound = feedbackSound;
end

%initialize event vector
bci=bciGenerateEvents(bci,'init');

% memory allocation
if isfield(bci.paradigm,'stimVelocityToDestination')&& bci.paradigm.stimVelocityToDestination>0,
    timeForFeedbackMovement = 10/bci.paradigm.stimVelocityToDestination*...
        sum(ismember(bci.events,[bci.leftEvent,bci.rightEvent,bci.baseEvent]));
else
    timeForFeedbackMovement=0;
end
trainRunDur = sum(bci.eventDur) + timeForFeedbackMovement + bci.preRunDur;
nData = round((max(bci.testRunDur,trainRunDur)+10) * bci.srate);
nEpochs = round((max(bci.testRunDur,trainRunDur)-length(bci.events)*bci.RT+10)*bci.init.epochRate);

nPreprocChans = length(bci.preprocChan);
if (size(MEMdata,1)~=bci.nChan || size(MEMdata,2)~=nData),
    MEMdata = zeros(bci.nChan,nData,'single');
end
if size(MEMepoch,1)~=nPreprocChans || size(MEMepoch,2)~=bci.init.nFeatPerChan || size(MEMepoch,3)~=nEpochs,
    MEMepoch = zeros(nPreprocChans,bci.init.nFeatPerChan,nEpochs);
end
